
COMPILE_OPT IDL2, LOGICAL_PREDICATE
COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile

conf = omni_load_conf()
s = omni_read_cat(conf.cat)
v_std = findgen(conf.nvbin)*conf.deltav + conf.vstart

IF n_elements(grs) EQ 0 THEN $
   restore,'./local/'+conf.survey+'_grs_spectra_r46.sav',/ver
restore,'./local/'+conf.survey+'_velocities.sav',/ver


;; Find where GRS spectrum exists
ig = where(total(grs.onspec,1) NE 0.,ng)
print,ng

;; Find where HCO+ or NH3 velocities exist
ih = where(v.mol[0].nspec NE 0, nh)
in = where(v.mol[3].nspec NE 0, nn)

print,nh,nn




;;=========================================================================
;; Plots


;;***********************************************
;; V-V plots
myps,'./grs_paper/brightco_comp.eps',xsize=3.1,ysize=1.75,thick=1.5

print,!x.omargin,!y.omargin
!x.omargin = [0,-2.5]
!y.omargin = [+0.25,-1.0]
multiplot,[2,1]

csf = 0.55
syms = 0.2

;;==========================
;; Top Left Panel

;; Get list of GRS-HCOP compatabilities
iginh = where_array(ig,ih,ngh)
ihing = where_array(ih,ig,nhg)
print,ngh,nhg
ihp = ih[iginh]
igp = ig[ihing]

;; Find vlsr of peak of ONSPEC
onspec  = grs[igp].onspec            ; ONSPEC (2-D)
v_arr   = cmreplicate(v_std,ngh)     ; Corresponding 2-D V_STD array
maxvals = max(onspec,jj,DIMENSION=1) ; Locations of maxima
gvlsr   = v_arr[jj]                  ; VLSR from GRS Brightest Peak
;; Set VLSR = -1000 if maxval < 0.1 K.
print,'MAXVALS: ',m4_stat(maxvals)
noni = where(maxvals LT 0.1, nnon)
print,'N below T_A = 0.1 K: ',nnon
IF nnon NE 0 THEN gvlsr[noni] = -1000.

YRANGE = [-10,145]

cgPlot,xr=[-50,130],yr=YRANGE,v[ihp].mol[0].vlsr,gvlsr,psym=16,$
       symsize=syms,/xst,$
       ytit='Brightest !u13!nCO V!dLSR!n  [km s!u-1!n]',charsize=1.0*csf,/yst,$
       xtit='HCO!u+!n  V!dLSR!n  [km s!u-1!n]',ytickinterval=50,/nodata

;; Draw lines in the background
cgColorfill,[!y.crange[0]+5,!x.crange[1],!x.crange[1],!y.crange[0]-5],$
            [!y.crange[0],!x.crange[1]-5,!x.crange[1]+5,!y.crange[0]],$
            color='sky blue'
vline,-5,color='blk5'
        

;; Plot black points
ipos = where(v[ihp].mol[0].vlsr GE -5, npos)
cgOplot,v[ihp[ipos]].mol[0].vlsr,gvlsr[ipos],$
        psym=16,symsize=syms,color='black'

;; Plot pink points
ineg = where(v[ihp].mol[0].vlsr LT -5, nneg)
cgOplot,v[ihp[ineg]].mol[0].vlsr,gvlsr[ineg],$
        psym=16,symsize=syms,color='Magenta'


;; Compute matching to w/in 5 km/s
ikms = where(abs(v[ihp].mol[0].vlsr-gvlsr) LE 5., nkms)
mfr_hp = float(nkms)/n_elements(ihp)*100.
al_legend,/top,/left,/clear,box=0,charsize=0.9*csf,$
          [string(mfr_hp,format="(F0.1)")+'%']


delv_hg = abs(v[ihp].mol[0].vlsr-gvlsr)
print,'LOOK HERE -- LOOK HERE -- LOOK HERE -- LOOK HERE -- LOOK HERE --'
print,m4_stat(double(v[ihp].mol[0].vlsr))
print,m4_stat(double(gvlsr))
print,m4_stat(double(delv_hg))

al_legend,/bottom,/right,box=0,['N = '+string(nhg,format="(I0)")],$
          charsize=0.9*csf

;; Redraw axes for pretty
cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,/yst,ytickformat='blank_axis',ytickinterval=50
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',ytickinterval=50


multiplot,/doxaxis

;;==========================
;; Top Right Panel


;; Get list of GRS-NH3 compatabilities
iginn = where_array(ig,in,ngn)
ining = where_array(in,ig,nng)
print,ngn,nng
inp = in[iginn]
igp = ig[ining]

;; Find vlsr of peak of ONSPEC
onspec  = grs[igp].onspec            ; ONSPEC (2-D)
v_arr   = cmreplicate(v_std,ngn)     ; Corresponding 2-D V_STD array
maxvals = max(onspec,jj,DIMENSION=1) ; Locations of maxima
gvlsr   = v_arr[jj]                  ; VLSR from GRS Brightest Peak

print,m4_stat(v[inp].cnum)
print,m4_stat(grs[igp].cnum)

cgPlot,xr=[-50,130],yr=YRANGE,v[inp].mol[3].vlsr,gvlsr,psym=16,$
       symsize=syms,/xst,charsize=1.0*csf,/yst,$
       xtit='NH!d3!n  V!dLSR!n  [km s!u-1!n]',ytickinterval=50,/nodata

;; Draw lines in the background
cgColorfill,[!y.crange[0]+5,!x.crange[1],!x.crange[1],!y.crange[0]-5],$
            [!y.crange[0],!x.crange[1]-5,!x.crange[1]+5,!y.crange[0]],$
            color='sky blue'
vline,-5,color='blk5'

;; Plot black points
cgOplot,v[inp].mol[3].vlsr,gvlsr,psym=16,symsize=syms

;; Compute matching to w/in 5 km/s
ikms = where(abs(v[inp].mol[3].vlsr-gvlsr) LE 5., nkms)
al_legend,/top,/left,/clear,box=0,charsize=0.9*csf,$
          [string(float(nkms)/n_elements(inp)*100.,format="(F0.1)")+'%']

al_legend,/bottom,/right,box=0,['N = '+string(nng,format="(I0)")],$
          charsize=0.9*csf

;; Redraw axes for pretty
cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,/yst,ytickformat='blank_axis',ytickinterval=50
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',ytickinterval=50


myps,/done,/mp


;;============================================
;; HCO+ - GRS Velocity comparison vs. S_{1.1}


myps,'./grs_paper/brightco_flux.eps',xsize=10,ysize=4.5
multiplot_xm,[2,1],mpcharsize=1.0,xgap=0.035

cgPlot,s[ihp].flux40,delv_hg,psym=16,symsize=syms,xtit='S!d40!n  [Jy]',$
       ytit='|'+cgSymbol('Delta')+'V!dLSR!n|  [km s!u-1!n]',$
       charsize=1.0,/xlog;,/ylog,yr=[1.d-6,1.d4]
vline,/h,5,/log,color='green'

al_legend,/top,/right,box=0,['N = '+string(n_elements(ihp),format="(I0)")]


multiplot,/doyaxis

;;======
;; Matching % as function of S40

yarr = histogram(alog10(s[ihp].flux40),bin=0.1,min=-2.0,max = 2.0,loc=xarr,$
                 reverse_indices=rind) ; RIND indexes IHP!
nx = n_elements(xarr)
dx = arrstep(xarr)
mfr = dblarr(nx)

FOR ii=0,nx-1 DO BEGIN
   IF rind[ii] EQ rind[ii+1] THEN CONTINUE
   
   inds    = rind[rind[ii] : rind[ii+1]-1]
   ind_ihp = ihp[inds]
   
   ;; Find fraction of delv_hg that are above 5 in each bin!
   ibad = where(delv_hg[inds] GT 5., nbad)
   ntot = n_elements(inds)
   
   mfr[ii] = double(ntot-nbad)/double(ntot)
   ;; print,xarr[ii],double(ntot-nbad)/double(ntot)
      
ENDFOR

cgPlot,xarr+dx,mfr,psym=10,charsize=1.0,yr=[0,1.05],xst=4,$
       ytit='Velocity Matching Fraction'
vline,/h,mfr_hp/100.,color='cyan'

cgAxis,xaxis=0,xr=10.^!x.crange,charsize=1.0,/save,xtit='S!d40!n  [Jy]',/xlog
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'

myps,/done,/mp














;;==========================
;; Bottom Left Panel
myps,'./grs_paper/mestco_comp.eps',xsize=3.1,ysize=1.75,thick=1.5

print,!x.omargin,!y.omargin
multiplot,[2,1]

print,m4_stat(v.grs.flag)

;; Update ihp to remove grs.flag = 0
igz = where(v[ihp].grs.flag NE 0 AND v[ihp].grs.flag NE 3, ngz)
print,ngz
ihp = ihp[igz]

cgPlot,xr=[-50,130],yr=YRANGE,v[ihp].vlsr,v[ihp].grs.vlsr,psym=16,$
       symsize=syms,/xst,/nodata,$
       ytit='Morph. Extracted !u13!nCO V!dLSR!n  [km s!u-1!n]',$
       charsize=1.0*csf, xtit='HCO!u+!n  V!dLSR!n  [km s!u-1!n]',$
       ytickinterval=50

;; Draw lines in the background
cgColorfill,[!y.crange[0]+5,!x.crange[1],!x.crange[1],!y.crange[0]-5],$
            [!y.crange[0],!x.crange[1]-5,!x.crange[1]+5,!y.crange[0]],$
            color='sky blue'
vline,-5,color='blk5'
        

;; Plot black points
ipos = where(v[ihp].mol[0].vlsr GE -5, npos)
cgOplot,v[ihp[ipos]].vlsr,v[ihp[ipos]].grs.vlsr,$
        psym=16,symsize=syms,color='black'

;; Plot pink points
ineg = where(v[ihp].mol[0].vlsr LT -5, nneg)
cgOplot,v[ihp[ineg]].vlsr,v[ihp[ineg]].grs.vlsr,$
        psym=16,symsize=syms,color='Magenta'

;; Compute matching to w/in 5 km/s
ikms = where(abs(v[ihp].mol[0].vlsr-v[ihp].grs.vlsr) LE 5., nkms)
ineg = where(v[ihp].mol[0].vlsr LE -5., nneg)
print,'DESIRED STATS (1/25/14):'
print,n_elements(ihp),nkms,n_elements(ihp)-nkms,nneg
mfr_hp = float(nkms)/n_elements(ihp)*100.
al_legend,/top,/left,/clear,box=0,charsize=0.9*csf,$
          [string(float(nkms)/n_elements(ihp)*100.,format="(F0.1)")+'%']
delv_hg = abs(v[ihp].mol[0].vlsr-v[ihp].grs.vlsr)

print,'LOOK HERE -- LOOK HERE -- LOOK HERE -- LOOK HERE -- LOOK HERE --'
print,m4_stat(double(v[ihp].mol[0].vlsr))
print,m4_stat(double(v[ihp].grs.vlsr))
print,m4_stat(double(delv_hg))

al_legend,/bottom,/right,box=0,['N = '+string(ngz,format="(I0)")],$
          charsize=0.9*csf

;; Redraw axes for pretty
cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,/yst,ytickformat='blank_axis',ytickinterval=50
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',ytickinterval=50


multiplot

;;==========================
;; Bottom Right Panel

;; Update inp to remove grs.flag = 0
igz = where(v[inp].grs.flag NE 0 AND v[inp].grs.flag NE 3, ngz)
print,ngz
inp = inp[igz]

cgPlot,xr=[-50,130],yr=YRANGE,v[inp].vlsr,v[inp].grs.vlsr,psym=16,$
       symsize=syms,/xst,/nodata,$
       xtit='NH!d3!n  V!dLSR!n  [km s!u-1!n]',charsize=1.0*csf,ytickinterval=50

;; Draw lines in the background
cgColorfill,[!y.crange[0]+5,!x.crange[1],!x.crange[1],!y.crange[0]-5],$
            [!y.crange[0],!x.crange[1]-5,!x.crange[1]+5,!y.crange[0]],$
            color='sky blue'
vline,-5,color='blk5'

;; Plot black points
cgOplot,v[inp].vlsr,v[inp].grs.vlsr,psym=16,symsize=syms

;; Compute matching to w/in 5 km/s
ikms = where(abs(v[inp].mol[3].vlsr-v[inp].grs.vlsr) LE 5., nkms)
al_legend,/top,/left,/clear,box=0,charsize=0.9*csf,$
          [string(float(nkms)/n_elements(inp)*100.,format="(F0.1)")+'%']

al_legend,/bottom,/right,box=0,['N = '+string(ngz,format="(I0)")],$
          charsize=0.9*csf

;; Redraw axes for pretty
cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
cgAxis,yaxis=0,/yst,ytickformat='blank_axis',ytickinterval=50
cgAxis,yaxis=1,/yst,ytickformat='blank_axis',ytickinterval=50

myps,/done,/mp



;;============================================
;; HCO+ - GRS Velocity comparison vs. S_{1.1}


myps,'./grs_paper/mestco_flux.eps',xsize=10,ysize=4.5
multiplot_xm,[2,1],mpcharsize=1.0,xgap=0.035

cgPlot,s[ihp].flux40,delv_hg,psym=16,symsize=syms,xtit='S!d40!n  [Jy]',$
       ytit='|'+cgSymbol('Delta')+'V!dLSR!n|  [km s!u-1!n]',$
       charsize=1.0,/xlog;,/ylog
al_legend,/top,/right,box=0,['N = '+string(n_elements(ihp),format="(I0)")]
vline,/h,5,/log,color='green'
multiplot,/doyaxis

;;======
;; Matching % as function of S40

yarr = histogram(alog10(s[ihp].flux40),bin=0.1,min=-2.0,max = 2.0,loc=xarr,$
                 reverse_indices=rind) ; RIND indexes IHP!
nx = n_elements(xarr)
dx = arrstep(xarr)
mfr = dblarr(nx)

FOR ii=0,nx-1 DO BEGIN
   IF rind[ii] EQ rind[ii+1] THEN CONTINUE
   
   inds    = rind[rind[ii] : rind[ii+1]-1]
   ind_ihp = ihp[inds]
   
   ;; Find fraction of delv_hg that are above 5 in each bin!
   ibad = where(delv_hg[inds] GT 5., nbad)
   ntot = n_elements(inds)
   
   mfr[ii] = double(ntot-nbad)/double(ntot)
   ;; print,xarr[ii],double(ntot-nbad)/double(ntot)
      
ENDFOR

cgPlot,xarr+dx,mfr,psym=10,charsize=1.0,yr=[0,1.05],xst=4,$
       ytit='Velocity Matching Fraction'
vline,/h,mfr_hp/100.,color='cyan'

cgAxis,xaxis=0,xr=10.^!x.crange,charsize=1.0,/save,xtit='S!d40!n  [Jy]',/xlog
cgAxis,xaxis=1,/xst,xtickformat='blank_axis'

myps,/done,/mp


;; ;;***********************************************
;; ;;***********************************************
;; ;;***********************************************
;; ;;***********************************************
;; ;;***********************************************
;; ;; Delta-V Histograms
;; myps,'./grs_paper/brightco_histo.eps',xsize=10,ysize=9

;; multiplot_xm,[2,2],mpcharsize=1.0,/rowmajor,xgap=0.02,$
;;              mytitle='N per bin'

;; ;;==========================
;; ;; Top Left Panel

;; ;; Get list of GRS-HCOP compatabilities
;; iginh = where_array(ig,ih,ngh)
;; ihing = where_array(ih,ig,nhg)
;; print,ngh,nhg
;; ihp = ih[iginh]
;; igp = ig[ihing]

;; ;; Find vlsr of peak of ONSPEC
;; onspec  = grs[igp].onspec            ; ONSPEC (2-D)
;; v_arr   = cmreplicate(v_std,ngh)     ; Corresponding 2-D V_STD array
;; maxvals = max(onspec,jj,DIMENSION=1) ; Locations of maxima
;; gvlsr   = v_arr[jj]                  ; VLSR from GRS Brightest Peak

;; plothist,v[ihp].mol[0].vlsr-gvlsr,charsize=1.0,xr=[-12,12],bin=0.2,$
;;          ytit='Brightest !u`3!nCO',xarr,yarr

;; yf = mpfitpeak(xarr,yarr,A,nterms=3)
;; cgOplot,xarr,yf,color='cyan'

;; al_legend,/top,/right,charsize=0.9,linestyle=0,color='cyan',linsize=0.4,$
;;           [cgSymbol('sigma')+' = '+string(A[2],format="(F0.2)")+' km s!u-1!n']
 


;; al_legend,/bottom,/right,box=0,['N = '+string(nhg,format="(I0)")]


;; multiplot,/doxaxis

;; ;;==========================
;; ;; Bottom Left Panel

;; ;; Update ihp to remove grs.flag = 0
;; igz = where(v[ihp].grs.flag NE 0 AND v[ihp].grs.flag NE 3, ngz)
;; print,ngz
;; ihp = ihp[igz]

;; plothist,v[ihp].vlsr-v[ihp].grs.vlsr,charsize=1.0,xr=[-12,12],bin=0.2,$
;;          ytit='On-Off !u`3!nCO',xarr,yarr,$
;;          xtit='HCO!u+!n V!dLSR!n - !u13!nCO V!dLSR!n  [km s!u-1!n]'

;; yf = mpfitpeak(xarr,yarr,A,nterms=3)
;; cgOplot,xarr,yf,color='cyan'

;; al_legend,/top,/right,charsize=0.9,linestyle=0,color='cyan',linsize=0.4,$
;;           [cgSymbol('sigma')+' = '+string(A[2],format="(F0.2)")+' km s!u-1!n']
 


;; al_legend,/bottom,/right,box=0,['N = '+string(ngz,format="(I0)")]

;; cgAxis,xaxis=0,/xst,xtickformat='blank_axis'

;; multiplot,/doyaxis

;; ;;==========================
;; ;; Top Right Panel


;; ;; Get list of GRS-NH3 compatabilities
;; iginn = where_array(ig,in,ngn)
;; ining = where_array(in,ig,nng)
;; print,ngn,nng
;; inp = in[iginn]
;; igp = ig[ining]

;; ;; Find vlsr of peak of ONSPEC
;; onspec  = grs[igp].onspec            ; ONSPEC (2-D)
;; v_arr   = cmreplicate(v_std,ngn)     ; Corresponding 2-D V_STD array
;; maxvals = max(onspec,jj,DIMENSION=1) ; Locations of maxima
;; gvlsr   = v_arr[jj]                  ; VLSR from GRS Brightest Peak

;; print,m4_stat(v[inp].cnum)
;; print,m4_stat(grs[igp].cnum)

;; plothist,v[inp].mol[3].vlsr-gvlsr,charsize=1.0,xr=[-12,12],bin=0.2,xarr,yarr


;; yf = mpfitpeak(xarr,yarr,A,nterms=3)
;; cgOplot,xarr,yf,color='cyan'

;; al_legend,/top,/right,charsize=0.9,linestyle=0,color='cyan',linsize=0.4,$
;;           [cgSymbol('sigma')+' = '+string(A[2],format="(F0.2)")+' km s!u-1!n']
 

;; al_legend,/bottom,/right,box=0,['N = '+string(nng,format="(I0)")]


;; multiplot,/doxaxis,/doyaxis

;; ;;==========================
;; ;; Bottom Right Panel

;; ;; Update inp to remove grs.flag = 0
;; igz = where(v[inp].grs.flag NE 0 AND v[inp].grs.flag NE 3, ngz)
;; print,ngz
;; inp = inp[igz]

;; plothist,v[inp].vlsr-v[inp].grs.vlsr,charsize=1.0,xr=[-12,12],bin=0.2,xarr,$
;;          xtit='NH!d3!n V!dLSR!n - !u13!nCO V!dLSR!n  [km s!u-1!n]',yarr


;; yf = mpfitpeak(xarr,yarr,A,nterms=3)
;; cgOplot,xarr,yf,color='cyan'

;; al_legend,/top,/right,charsize=0.9,linestyle=0,color='cyan',linsize=0.4,$
;;           [cgSymbol('sigma')+' = '+string(A[2],format="(F0.2)")+' km s!u-1!n']
 



;; al_legend,/bottom,/right,box=0,['N = '+string(ngz,format="(I0)")]

;; cgAxis,xaxis=0,/xst,xtickformat='blank_axis'


;; myps,/done,/mp
END












